from libmyatomism import *
#from math import *
#import array 
#import matplotlib.pyplot as plt

msLogger.setPriorityMin(Priority.ERROR)

unit=msUnitsManager.New("Angstrom Degree amu kcal/mol ps")

q0 = msScalarVariable.New("Angstrom").set(5,0.5,6,0.1,0.1).setId("q0")
q1 = msScalarVariable.New("Degree").set(120,0,360,2,2).setId("q1")
q2 = msScalarVariable.New("Degree").set(120,0,360,2,2).setId("q2")
q3 = msScalarVariable.New("Degree").set(120,0,360,2,2).setId("q3")

gcoors = msGeneralizedCoordinates.New(unit).addVar(q0).addVar(q1).addVar(q2).addVar(q3)

zmat = msZMat.New(unit).setId("zmat").set("N N 0 1.1 N 1 5. 0 10. N 2 1.1 1 90. 0 0.")

kinfct0 = msParser.New("Angstrom").setCoordinates(gcoors).setExpression("q0").setId("f=q0")
kinfct1 = msParser.New("Degree").setCoordinates(gcoors).setExpression("q1").setId("f=q1")
kinfct2 = msParser.New("Degree").setCoordinates(gcoors).setExpression("q2").setId("f=q2")
kinfct3 = msParser.New("Degree").setCoordinates(gcoors).setExpression("q3").setId("f=q3")

system=msSystem.New(unit).addEntity(zmat)

kinop=msKineticOperator.New(unit).set(gcoors,system).addDynamicDof(zmat,1,kinfct0).addDynamicDof(zmat,2,kinfct1).addDynamicDof(zmat,4,kinfct2).addDynamicDof(zmat,5,kinfct3)

cfour = qcCfour.New().setLevel( BasisSet.ANO0 , CiType.CCSDpT ).setGenBasSymLink("/org/centers/pecos/APPLICATIONS/cfour/basis/GENBAS")

zmatN2 = msZMat.New(unit).setId("zmat").set("N N 0 1.1")
eN2    = cfour.singlePoint(zmatN2)

#view=msVisualizerAtoms.New()
#view.watch(system)
file = open('approach0.txt','w')
gcoors.setCoordinate(1,120)
gcoors.setCoordinate(2,120)
gcoors.setCoordinate(3,120.)
for i in range(0,10):
  gcoors.setCoordinate(0,6-4*i/10.)
  kinop.setDynamicDofs()
  e = cfour.singlePoint(zmat)
  print q0.Value,"\t",q1.Value,"\t",q2.Value,"\t",q3.Value,"\t",e,"\t",unit.convert("hartree",e-2*eN2)
  file.write(str(q0.Value)+" "+str(q1.Value)+" "+str(q2.Value)+" "+str(q3.Value)+" "+str(e)+" "+str(unit.convert("hartree",e-2*eN2))+"\n" )
file.close()
  
file = open('datafile.txt','w')

for i in range(0,1000):
  sepmin = 0
  while( sepmin < 1 ):
    gcoors.setRandom()
    kinop.setDynamicDofs()
    sepmin = system.minSeparation()
    e = cfour.singlePoint(zmat)
  print q0.Value,"\t",q1.Value,"\t",q2.Value,"\t",q3.Value,"\t",e,"\t",unit.convert("hartree",e-2*eN2)
  file.write(str(q0.Value)+" "+str(q1.Value)+" "+str(q2.Value)+" "+str(q3.Value)+" "+str(e)+" "+str(unit.convert("hartree",e-2*eN2))+"\n" )

file.close()